{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# System building: Protein in Membrane with Ligand"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "In this tutorial, we will showcase how to build a protein ligand system for simulating binding. The sample system is Trypsin (the protein) and benzamidine (the ligand).\n",
    "\n",
    "Let's start by doing some imports and definitions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Please cite HTMD: Doerr et al.(2016)JCTC,12,1845. \n",
      "https://dx.doi.org/10.1021/acs.jctc.6b00049\n",
      "Documentation: http://software.acellera.com/\n",
      "To update: conda update htmd -c acellera -c psi4\n",
      "\n",
      "You are on the latest HTMD version (unpackaged : /home/joao/maindisk/software/repos/cuzzo87/htmd/htmd).\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from htmd.ui import *\n",
    "from htmd.home import home\n",
    "from os.path import join\n",
    "config(viewer='webgl')\n",
    "datadir = home(dataDir='building-protein-ligand')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Load the protein-ligand complex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can obtain the protein-ligand complex from the PDB database (ID:3PTB). The complex is already available in the data distributed with HTMD and either one could be used:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-04-30 14:51:38,239 - htmd.molecule.readers - INFO - Using local copy for 3PTB: /home/joao/maindisk/software/repos/cuzzo87/htmd/htmd/data/pdb/3ptb.pdb\n",
      "2018-04-30 14:51:38,361 - htmd.molecule.molecule - WARNING - Residue insertions were detected in the Molecule. It is recommended to renumber the residues using the Molecule.renumberResidues() method.\n",
      "2018-04-30 14:51:38,473 - htmd.molecule.molecule - WARNING - Residue insertions were detected in the Molecule. It is recommended to renumber the residues using the Molecule.renumberResidues() method.\n"
     ]
    }
   ],
   "source": [
    "# One can download it directly from the RCSB servers\n",
    "prot = Molecule('3PTB')\n",
    "# Or use the pdb file found in the HTMD data directory\n",
    "prot = Molecule(join(datadir, 'trypsin.pdb'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "326c6f1c1fd140408f489101f575f0d6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "prot.view()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Clean the structures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "The PDB crystal structure contains the protein as well as water molecules, a calcium ion and a ligand. Here we will start by removing the ligand from the protein Molecule as we will add it later to manipulate it separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-04-30 14:51:38,911 - htmd.molecule.molecule - INFO - Removed 9 atoms. 1692 atoms remaining in the molecule.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([1630, 1631, 1632, 1633, 1634, 1635, 1636, 1637, 1638], dtype=int32)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prot.remove('resname BEN')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Preparing the protein"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this step, we prepare the protein for simulation by adding hydrogens, setting the protonation states, and optimizing the protein (more details on the protein preparation tutorial):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-04-30 14:51:39,057 - propka - INFO - No pdbfile provided\n",
      "2018-04-30 14:51:41,578 - htmd.builder.preparation - WARNING - The following residue has not been optimized: CA\n",
      "2018-04-30 14:51:49,983 - htmd.builder.preparationdata - INFO - The following residues are in a non-standard state: CYS     6  A (CYX), HIS    22  A (HIE), CYS    24  A (CYX), HIS    39  A (HIP), CYS    40  A (CYX), HIS    72  A (HID), CYS   108  A (CYX), CYS   115  A (CYX), CYS   136  A (CYX), CYS   147  A (CYX), CYS   161  A (CYX), CYS   172  A (CYX), CYS   182  A (CYX), CYS   196  A (CYX), CYS   209  A (CYX)\n",
      "2018-04-30 14:51:50,022 - htmd.builder.preparationdata - WARNING - Dubious protonation state: the pKa of 4 residues is within 1.0 units of pH 7.0.\n",
      "2018-04-30 14:51:50,026 - htmd.builder.preparationdata - WARNING - Dubious protonation state:    HIS    39  A (pKa= 7.46)\n",
      "2018-04-30 14:51:50,029 - htmd.builder.preparationdata - WARNING - Dubious protonation state:    GLU    51  A (pKa= 6.10)\n",
      "2018-04-30 14:51:50,030 - htmd.builder.preparationdata - WARNING - Dubious protonation state:    ASP   170  A (pKa= 6.49)\n",
      "2018-04-30 14:51:50,031 - htmd.builder.preparationdata - WARNING - Dubious protonation state:    N+      0T A (pKa= 7.49)\n",
      "2018-04-30 14:51:50,121 - htmd.builder.preparationdata - WARNING - Found N-terminus 80.7% buried (> 50.0% threshold)\n",
      "2018-04-30 14:51:50,122 - htmd.builder.preparationdata - WARNING - Found C-terminus involved in H bonds\n"
     ]
    }
   ],
   "source": [
    "prot = proteinPrepare(prot, pH=7.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Define segments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "To build a system in HTMD, we need to separate the chemical molecules into separate segments. This prevents the builder from accidentally bonding different chemical molecules and allows us to add caps to them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "prot = autoSegment(prot, sel='protein')\n",
    "prot.set('segid', 'W', sel='water')\n",
    "prot.set('segid', 'CA', sel='resname CA')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Center the protein to the origin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "prot.center()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Let's work on the ligand!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the ligand from the HTMD data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-04-30 14:51:50,585 - htmd.molecule.readers - WARNING - Element of atom ID 1 could not be automatically guessed from its MOL2 atomtype (ca).\n",
      "2018-04-30 14:51:50,586 - htmd.molecule.readers - WARNING - Element of atom ID 2 could not be automatically guessed from its MOL2 atomtype (ca).\n",
      "2018-04-30 14:51:50,587 - htmd.molecule.readers - WARNING - Element of atom ID 3 could not be automatically guessed from its MOL2 atomtype (ca).\n",
      "2018-04-30 14:51:50,588 - htmd.molecule.readers - WARNING - Element of atom ID 4 could not be automatically guessed from its MOL2 atomtype (ca).\n",
      "2018-04-30 14:51:50,589 - htmd.molecule.readers - WARNING - Element of atom ID 5 could not be automatically guessed from its MOL2 atomtype (ca).\n",
      "2018-04-30 14:51:50,590 - htmd.molecule.readers - WARNING - Element of atom ID 6 could not be automatically guessed from its MOL2 atomtype (ca).\n",
      "2018-04-30 14:51:50,592 - htmd.molecule.readers - WARNING - Element of atom ID 7 could not be automatically guessed from its MOL2 atomtype (ce).\n",
      "2018-04-30 14:51:50,593 - htmd.molecule.readers - WARNING - Element of atom ID 8 could not be automatically guessed from its MOL2 atomtype (ha).\n",
      "2018-04-30 14:51:50,594 - htmd.molecule.readers - WARNING - Element of atom ID 9 could not be automatically guessed from its MOL2 atomtype (ha).\n",
      "2018-04-30 14:51:50,595 - htmd.molecule.readers - WARNING - Element of atom ID 10 could not be automatically guessed from its MOL2 atomtype (ha).\n",
      "2018-04-30 14:51:50,596 - htmd.molecule.readers - WARNING - Element of atom ID 11 could not be automatically guessed from its MOL2 atomtype (ha).\n",
      "2018-04-30 14:51:50,597 - htmd.molecule.readers - WARNING - Element of atom ID 12 could not be automatically guessed from its MOL2 atomtype (ha).\n",
      "2018-04-30 14:51:50,598 - htmd.molecule.readers - WARNING - Element of atom ID 13 could not be automatically guessed from its MOL2 atomtype (nh).\n",
      "2018-04-30 14:51:50,599 - htmd.molecule.readers - WARNING - Element of atom ID 14 could not be automatically guessed from its MOL2 atomtype (nh).\n",
      "2018-04-30 14:51:50,600 - htmd.molecule.readers - WARNING - Element of atom ID 15 could not be automatically guessed from its MOL2 atomtype (hn).\n",
      "2018-04-30 14:51:50,601 - htmd.molecule.readers - WARNING - Element of atom ID 16 could not be automatically guessed from its MOL2 atomtype (hn).\n",
      "2018-04-30 14:51:50,602 - htmd.molecule.readers - WARNING - Element of atom ID 17 could not be automatically guessed from its MOL2 atomtype (hn).\n",
      "2018-04-30 14:51:50,603 - htmd.molecule.readers - WARNING - Element of atom ID 18 could not be automatically guessed from its MOL2 atomtype (hn).\n"
     ]
    }
   ],
   "source": [
    "ligand = Molecule(join(datadir, 'benzamidine.mol2'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's center the ligand and visualize it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "290df5adc31f4790af97103b1223b5c4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ligand.center()\n",
    "ligand.view()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can give a convenient segid and resname to the ligand\n",
    "# The resname should be MOL to match the parameters in the\n",
    "# rtf and prm files.\n",
    "ligand.set('segid','L')\n",
    "ligand.set('resname','MOL')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But the ligand is now located inside the protein...\n",
    "We would like the ligand to be:\n",
    "\n",
    "* At a certain distance from the protein\n",
    "* Rotated randomly, to provide different starting conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Let's randomize the ligand position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "ligand.rotateBy(uniformRandomRotation())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This took care of the ligand rotation around its own center. \n",
    "We still need to position it far from the protein.\n",
    "First, find out the radius of the protein:\n",
    "\n",
    "![maxdist](http://pub.htmd.org/tutorials/system-building-protein-ligand/maxdist.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "28.832803412\n"
     ]
    }
   ],
   "source": [
    "from moleculekit.util import maxDistance\n",
    "D = maxDistance(prot, 'all')\n",
    "print(D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "D += 10\n",
    "# Move the ligand 10 Angstrom away from the furthest protein atom in X dimension\n",
    "ligand.moveBy([D, 0, 0])  \n",
    "# rotateBy rotates by default around [0, 0, 0]. Since the ligand has been moved\n",
    "# away from the center it will be rotated in a sphere of radius D+10 around [0, 0, 0]\n",
    "ligand.rotateBy(uniformRandomRotation())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Mix it all together"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a00dbee7035f4390876f49f6b9c483bc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mol = Molecule(name='combo')\n",
    "mol.append(prot)\n",
    "mol.append(ligand)\n",
    "mol.reps.add(sel='protein', style='NewCartoon', color='Secondary Structure')\n",
    "mol.reps.add(sel='resname MOL', style='Licorice')\n",
    "mol.view()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Solvate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "> Water is the driving force of all nature. --Leonardo da Vinci\n",
    "\n",
    "![waterbox](http://pub.htmd.org/tutorials/system-building-protein-ligand/waterbox.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-04-30 14:51:51,230 - htmd.builder.solvate - INFO - Using water pdb file at: /data/joao/maindisk/software/repos/cuzzo87/htmd/htmd/builder/wat.pdb\n",
      "2018-04-30 14:51:52,328 - htmd.builder.solvate - INFO - Replicating 8 water segments, 2 by 2 by 2\n",
      "Solvating: 100%|██████████| 8/8 [00:09<00:00,  1.19s/it]\n",
      "2018-04-30 14:52:04,572 - htmd.builder.solvate - INFO - 20153 water molecules were added to the system.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "660916c6d30b4ba1b33655a89950570a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# We solvate with a larger box to fully solvate the ligand\n",
    "DW = D + 5\n",
    "smol = solvate(mol, minmax=[[-DW, -DW, -DW], [DW, DW, DW]])\n",
    "smol.reps.add(sel='water', style='Lines')\n",
    "smol.view()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Build the System with a specific forcefield"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "HTMD aims to be force-field agnostic. After you have a built system, you can either build it in Amber or CHARMM. The following sections work on the same previously solvated system and can be interconverted.\n",
    "\n",
    "Special care must be taken care in this case due to the use of benzamidine, which is not present by default on the respective forcefields."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### CHARMM forcefield"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---- Topologies files list: /home/joao/maindisk/software/repos/cuzzo87/htmd/htmd/builder/charmmfiles/top/ ----\n",
      "top/top_all22_prot.rtf\n",
      "top/top_all22star_prot.rtf\n",
      "top/top_all35_ethers.rtf\n",
      "top/top_all36_carb.rtf\n",
      "top/top_all36_cgenff.rtf\n",
      "top/top_all36_lipid.rtf\n",
      "top/top_all36_na.rtf\n",
      "top/top_all36_prot.rtf\n",
      "top/top_water_ions.rtf\n",
      "---- Parameters files list: /home/joao/maindisk/software/repos/cuzzo87/htmd/htmd/builder/charmmfiles/par/ ----\n",
      "par/par_all22_prot.prm\n",
      "par/par_all22star_prot.prm\n",
      "par/par_all35_ethers.prm\n",
      "par/par_all36_carb.prm\n",
      "par/par_all36_cgenff.prm\n",
      "par/par_all36_lipid.prm\n",
      "par/par_all36_na.prm\n",
      "par/par_all36_prot.prm\n",
      "par/par_all36_prot_mod.prm\n",
      "par/par_water_ions.prm\n",
      "---- Stream files list: /home/joao/maindisk/software/repos/cuzzo87/htmd/htmd/builder/charmmfiles/str/ ----\n",
      "str/carb/toppar_all36_carb_glycolipid.str\n",
      "str/carb/toppar_all36_carb_glycopeptide.str\n",
      "str/carb/toppar_all36_carb_model.str\n",
      "str/lipid/toppar_all36_lipid_bacterial.str\n",
      "str/lipid/toppar_all36_lipid_cardiolipin.str\n",
      "str/lipid/toppar_all36_lipid_cholesterol.str\n",
      "str/lipid/toppar_all36_lipid_cholesterol_model_1.str\n",
      "str/lipid/toppar_all36_lipid_detergent.str\n",
      "str/lipid/toppar_all36_lipid_inositol.str\n",
      "str/lipid/toppar_all36_lipid_list.str\n",
      "str/lipid/toppar_all36_lipid_miscellaneous.str\n",
      "str/lipid/toppar_all36_lipid_model.str\n",
      "str/lipid/toppar_all36_lipid_prot.str\n",
      "str/lipid/toppar_all36_lipid_sphingo.str\n",
      "str/lipid/toppar_all36_lipid_yeast.str\n",
      "str/misc/toppar_amines.str\n",
      "str/misc/toppar_dum_noble_gases.str\n",
      "str/misc/toppar_hbond.str\n",
      "str/misc/toppar_ions_won.str\n",
      "str/na/toppar_all36_na_model.str\n",
      "str/na/toppar_all36_na_modifications.str\n",
      "str/na/toppar_all36_na_nad_ppi.str\n",
      "str/na/toppar_all36_na_reactive_rna.str\n",
      "str/na/toppar_all36_na_rna_modified.str\n",
      "str/prot/toppar_all36_prot_aldehydes.str\n",
      "str/prot/toppar_all36_prot_arg0.str\n",
      "str/prot/toppar_all36_prot_d_aminoacids.str\n",
      "str/prot/toppar_all36_prot_fluoro_alkanes.str\n",
      "str/prot/toppar_all36_prot_heme.str\n",
      "str/prot/toppar_all36_prot_mod_d_aminoacids.str\n",
      "str/prot/toppar_all36_prot_model.str\n",
      "str/prot/toppar_all36_prot_modify_res.str\n",
      "str/prot/toppar_all36_prot_na_combined.str\n",
      "str/prot/toppar_all36_prot_pyridines.str\n",
      "str/prot/toppar_all36_prot_retinol.str\n"
     ]
    }
   ],
   "source": [
    "charmm.listFiles()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Build and ionize using CHARMM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-04-30 14:52:13,118 - htmd.builder.charmm - INFO - Writing out segments.\n",
      "2018-04-30 14:52:22,002 - htmd.builder.builder - INFO - 6 disulfide bonds were added\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bond between A: [serial 351 resid 24 resname CYX chain A segid P0]\n",
      "             B: [serial 571 resid 40 resname CYX chain A segid P0]\n",
      "\n",
      "Bond between A: [serial 1683 resid 115 resname CYX chain A segid P0]\n",
      "             B: [serial 2611 resid 182 resname CYX chain A segid P0]\n",
      "\n",
      "Bond between A: [serial 2146 resid 147 resname CYX chain A segid P0]\n",
      "             B: [serial 2353 resid 161 resname CYX chain A segid P0]\n",
      "\n",
      "Bond between A: [serial 1604 resid 108 resname CYX chain A segid P0]\n",
      "             B: [serial 3004 resid 209 resname CYX chain A segid P0]\n",
      "\n",
      "Bond between A: [serial 2494 resid 172 resname CYX chain A segid P0]\n",
      "             B: [serial 2799 resid 196 resname CYX chain A segid P0]\n",
      "\n",
      "Bond between A: [serial 92 resid 6 resname CYX chain A segid P0]\n",
      "             B: [serial 1988 resid 136 resname CYX chain A segid P0]\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-04-30 14:52:22,528 - htmd.builder.charmm - INFO - Starting the build.\n",
      "2018-04-30 14:52:23,368 - htmd.builder.charmm - WARNING - Failed to set coordinates for 8 atoms.\n",
      "2018-04-30 14:52:23,370 - htmd.builder.charmm - WARNING - Failed to guess coordinates for 3 atoms due to bad angles.\n",
      "2018-04-30 14:52:23,371 - htmd.builder.charmm - WARNING - Poorly guessed coordinates for 15 atoms.\n",
      "2018-04-30 14:52:23,372 - htmd.builder.charmm - WARNING - Please check /data/joao/maindisk/software/repos/cuzzo87/htmd/tutorials/build_charmm/log.txt for further information.\n",
      "2018-04-30 14:52:23,373 - htmd.builder.charmm - INFO - Finished building.\n",
      "2018-04-30 14:52:27,602 - htmd.builder.ionize - INFO - Adding 10 anions + 0 cations for neutralizing and 0 ions for the given salt concentration.\n",
      "2018-04-30 14:52:28,112 - htmd.builder.ionize - INFO - Min distance of ions from molecule: 5A\n",
      "2018-04-30 14:52:28,113 - htmd.builder.ionize - INFO - Min distance between ions: 5A\n",
      "2018-04-30 14:52:28,114 - htmd.builder.ionize - INFO - Placing 10 anions and 0 cations.\n",
      "2018-04-30 14:52:35,771 - htmd.builder.charmm - INFO - Writing out segments.\n",
      "2018-04-30 14:52:45,692 - htmd.builder.charmm - INFO - Starting the build.\n",
      "2018-04-30 14:52:46,518 - htmd.builder.charmm - WARNING - Failed to set coordinates for 8 atoms.\n",
      "2018-04-30 14:52:46,519 - htmd.builder.charmm - WARNING - Failed to guess coordinates for 3 atoms due to bad angles.\n",
      "2018-04-30 14:52:46,520 - htmd.builder.charmm - WARNING - Poorly guessed coordinates for 15 atoms.\n",
      "2018-04-30 14:52:46,521 - htmd.builder.charmm - WARNING - Please check /data/joao/maindisk/software/repos/cuzzo87/htmd/tutorials/build_charmm/log.txt for further information.\n",
      "2018-04-30 14:52:46,522 - htmd.builder.charmm - INFO - Finished building.\n"
     ]
    }
   ],
   "source": [
    "topos_charmm  = charmm.defaultTopo() + [join(datadir, 'benzamidine.rtf')]\n",
    "params_charmm = charmm.defaultParam() + [join(datadir, 'benzamidine.prm')]\n",
    "\n",
    "bmol_charmm = charmm.build(smol, topo=topos_charmm, param=params_charmm, outdir='./build_charmm')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### AMBER forcefield"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---- Forcefield files list: /home/joao/miniconda3/dat/leap/cmd/ ----\n",
      "leaprc.constph\n",
      "leaprc.DNA.bsc1\n",
      "leaprc.DNA.OL15\n",
      "leaprc.ff14SB.redq\n",
      "leaprc.ffAM1\n",
      "leaprc.ffPM3\n",
      "leaprc.gaff\n",
      "leaprc.gaff2\n",
      "leaprc.GLYCAM_06EPb\n",
      "leaprc.GLYCAM_06j-1\n",
      "leaprc.lipid14\n",
      "leaprc.modrna08\n",
      "leaprc.phosaa10\n",
      "leaprc.protein.fb15\n",
      "leaprc.protein.ff03.r1\n",
      "leaprc.protein.ff03ua\n",
      "leaprc.protein.ff14SB\n",
      "leaprc.protein.ff14SBonlysc\n",
      "leaprc.protein.ff15ipq\n",
      "leaprc.protein.ff15ipq-vac\n",
      "leaprc.RNA.OL3\n",
      "leaprc.RNA.YIL\n",
      "leaprc.water.opc\n",
      "leaprc.water.spce\n",
      "leaprc.water.spceb\n",
      "leaprc.water.tip3p\n",
      "leaprc.water.tip4pew\n",
      "leaprc.xFPchromophores\n",
      "---- OLD Forcefield files list: /home/joao/miniconda3/dat/leap/cmd/ ----\n",
      "oldff/leaprc.DNA.bsc0\n",
      "oldff/leaprc.ff02\n",
      "oldff/leaprc.ff02pol.r0\n",
      "oldff/leaprc.ff02pol.r1\n",
      "oldff/leaprc.ff02polEP.r0\n",
      "oldff/leaprc.ff02polEP.r1\n",
      "oldff/leaprc.ff03\n",
      "oldff/leaprc.ff10\n",
      "oldff/leaprc.ff14ipq\n",
      "oldff/leaprc.ff14SB\n",
      "oldff/leaprc.ff84\n",
      "oldff/leaprc.ff86\n",
      "oldff/leaprc.ff94\n",
      "oldff/leaprc.ff94.nmr\n",
      "oldff/leaprc.ff96\n",
      "oldff/leaprc.ff98\n",
      "oldff/leaprc.ff99\n",
      "oldff/leaprc.ff99bsc0\n",
      "oldff/leaprc.ff99SB\n",
      "oldff/leaprc.ff99SBildn\n",
      "oldff/leaprc.ff99SBnmr\n",
      "oldff/leaprc.GLYCAM_04\n",
      "oldff/leaprc.GLYCAM_06\n",
      "oldff/leaprc.GLYCAM_06EP\n",
      "oldff/leaprc.GLYCAM_06h\n",
      "oldff/leaprc.GLYCAM_06h-1\n",
      "oldff/leaprc.GLYCAM_06h-12SB\n",
      "oldff/leaprc.GLYCAM_06j_10\n",
      "oldff/leaprc.lipid11\n",
      "oldff/leaprc.parmbsc0_chiOL4_ezOL1\n",
      "oldff/leaprc.rna.ff02\n",
      "oldff/leaprc.rna.ff02EP\n",
      "oldff/leaprc.rna.ff84\n",
      "oldff/leaprc.rna.ff94\n",
      "oldff/leaprc.rna.ff98\n",
      "oldff/leaprc.rna.ff99\n",
      "oldff/leaprc.toyrna\n",
      "---- Topology files list: /home/joao/miniconda3/dat/leap/prep/ ----\n",
      "all_amino03.in\n",
      "all_aminoct03.in\n",
      "all_aminont03.in\n",
      "amino10.in\n",
      "amino12.in\n",
      "aminoct10.in\n",
      "aminoct12.in\n",
      "aminont10.in\n",
      "aminont12.in\n",
      "chcl3.in\n",
      "dna_nuc94-bsc0_chiOl4-ezOL1.in\n",
      "GLYCAM_06EPb.prep\n",
      "GLYCAM_06j-1.prep\n",
      "GLYCAM_lipids_06h.prep\n",
      "meoh.in\n",
      "nma.in\n",
      "nucleic10.in\n",
      "toyrna.in\n",
      "uni_amino03.in\n",
      "uni_aminoct03.in\n",
      "uni_aminont03.in\n",
      "---- Parameter files list: /home/joao/miniconda3/dat/leap/parm/ ----\n",
      "frcmod.chcl3\n",
      "frcmod.chiOL4\n",
      "frcmod.constph\n",
      "frcmod.dc4\n",
      "frcmod.DNA.OL15\n",
      "frcmod.fb15\n",
      "frcmod.ff02pol.r1\n",
      "frcmod.ff03\n",
      "frcmod.ff03ua\n",
      "frcmod.ff12SB\n",
      "frcmod.ff14SB\n",
      "frcmod.ff99bsc0CG\n",
      "frcmod.ff99SB\n",
      "frcmod.ff99SB14\n",
      "frcmod.ff99SBildn\n",
      "frcmod.ff99SBnmr\n",
      "frcmod.ff99SP\n",
      "frcmod.ions1lm_1264_spce\n",
      "frcmod.ions1lm_1264_tip3p\n",
      "frcmod.ions1lm_1264_tip4pew\n",
      "frcmod.ions1lm_126_hfe_opc\n",
      "frcmod.ions1lm_126_iod_opc\n",
      "frcmod.ions1lm_126_spce\n",
      "frcmod.ions1lm_126_tip3p\n",
      "frcmod.ions1lm_126_tip4pew\n",
      "frcmod.ions1lm_iod\n",
      "frcmod.ions234lm_1264_spce\n",
      "frcmod.ions234lm_1264_tip3p\n",
      "frcmod.ions234lm_1264_tip4pew\n",
      "frcmod.ions234lm_126_spce\n",
      "frcmod.ions234lm_126_tip3p\n",
      "frcmod.ions234lm_126_tip4pew\n",
      "frcmod.ions234lm_hfe_spce\n",
      "frcmod.ions234lm_hfe_tip3p\n",
      "frcmod.ions234lm_hfe_tip4pew\n",
      "frcmod.ions234lm_iod_spce\n",
      "frcmod.ions234lm_iod_tip3p\n",
      "frcmod.ions234lm_iod_tip4pew\n",
      "frcmod.ionsff99_tip3p\n",
      "frcmod.ionsjc_spce\n",
      "frcmod.ionsjc_tip3p\n",
      "frcmod.ionsjc_tip4pew\n",
      "frcmod.meoh\n",
      "frcmod.nma\n",
      "frcmod.opc\n",
      "frcmod.parmbsc0\n",
      "frcmod.parmbsc0_ezOL1\n",
      "frcmod.parmbsc1\n",
      "frcmod.parmCHI_YIL\n",
      "frcmod.phosaa10\n",
      "frcmod.pol3\n",
      "frcmod.protonated_nucleic\n",
      "frcmod.qspcfw\n",
      "frcmod.spce\n",
      "frcmod.spceb\n",
      "frcmod.spcfw\n",
      "frcmod.tip3pf\n",
      "frcmod.tip3pfb\n",
      "frcmod.tip4p\n",
      "frcmod.tip4pew\n",
      "frcmod.tip4pfb\n",
      "frcmod.tip5p\n",
      "frcmod.urea\n",
      "frcmod.vdWall\n",
      "frcmod.xFPchromophores\n",
      "---- Extra forcefield files list: /home/joao/maindisk/software/repos/cuzzo87/htmd/htmd/builder/amberfiles/ ----\n",
      "ff-nucleic-OL15/leaprc.ff-nucleic-OL15\n",
      "---- Extra topology files list: /home/joao/maindisk/software/repos/cuzzo87/htmd/htmd/builder/amberfiles/ ----\n",
      "cofactors/cofactors.in\n",
      "ff-ncaa/ffncaa.in\n",
      "ff-nucleic-OL15/leap-ff-nucleic-OL15.in\n",
      "ff-ptm/ffptm.in\n",
      "---- Extra parameter files list: /home/joao/maindisk/software/repos/cuzzo87/htmd/htmd/builder/amberfiles/ ----\n",
      "cofactors/ADP.frcmod\n",
      "cofactors/AMP.frcmod\n",
      "cofactors/ATP.frcmod\n",
      "cofactors/FMN.frcmod\n",
      "cofactors/GDP.frcmod\n",
      "cofactors/GTP.frcmod\n",
      "cofactors/HEM.frcmod\n",
      "cofactors/NAD.frcmod\n",
      "cofactors/NAI.frcmod\n",
      "cofactors/NAP.frcmod\n",
      "ff-ncaa/004.frcmod\n",
      "ff-ncaa/03Y.frcmod\n",
      "ff-ncaa/0A1.frcmod\n",
      "ff-ncaa/0BN.frcmod\n",
      "ff-ncaa/1MH.frcmod\n",
      "ff-ncaa/2AS.frcmod\n",
      "ff-ncaa/2GX.frcmod\n",
      "ff-ncaa/2ML.frcmod\n",
      "ff-ncaa/4IN.frcmod\n",
      "ff-ncaa/4PH.frcmod\n",
      "ff-ncaa/5JP.frcmod\n",
      "ff-ncaa/AA4.frcmod\n",
      "ff-ncaa/ABA.frcmod\n",
      "ff-ncaa/AHP.frcmod\n",
      "ff-ncaa/ALC.frcmod\n",
      "ff-ncaa/ALN.frcmod\n",
      "ff-ncaa/APD.frcmod\n",
      "ff-ncaa/BB8.frcmod\n",
      "ff-ncaa/BCS.frcmod\n",
      "ff-ncaa/CCS.frcmod\n",
      "ff-ncaa/CSA.frcmod\n",
      "ff-ncaa/D4P.frcmod\n",
      "ff-ncaa/DAB.frcmod\n",
      "ff-ncaa/DPP.frcmod\n",
      "ff-ncaa/ESC.frcmod\n",
      "ff-ncaa/FGL.frcmod\n",
      "ff-ncaa/GHG.frcmod\n",
      "ff-ncaa/GME.frcmod\n",
      "ff-ncaa/GNC.frcmod\n",
      "ff-ncaa/HHK.frcmod\n",
      "ff-ncaa/HLU.frcmod\n",
      "ff-ncaa/HLX.frcmod\n",
      "ff-ncaa/HOX.frcmod\n",
      "ff-ncaa/HPE.frcmod\n",
      "ff-ncaa/HQA.frcmod\n",
      "ff-ncaa/HTR.frcmod\n",
      "ff-ncaa/I2M.frcmod\n",
      "ff-ncaa/IGL.frcmod\n",
      "ff-ncaa/IIL.frcmod\n",
      "ff-ncaa/IML.frcmod\n",
      "ff-ncaa/KYN.frcmod\n",
      "ff-ncaa/LME.frcmod\n",
      "ff-ncaa/LMQ.frcmod\n",
      "ff-ncaa/ME0.frcmod\n",
      "ff-ncaa/MEA.frcmod\n",
      "ff-ncaa/MEN.frcmod\n",
      "ff-ncaa/MEQ.frcmod\n",
      "ff-ncaa/MLE.frcmod\n",
      "ff-ncaa/MLZ.frcmod\n",
      "ff-ncaa/MME.frcmod\n",
      "ff-ncaa/MMO.frcmod\n",
      "ff-ncaa/MVA.frcmod\n",
      "ff-ncaa/NAL.frcmod\n",
      "ff-ncaa/NCY.frcmod\n",
      "ff-ncaa/NLE.frcmod\n",
      "ff-ncaa/NVA.frcmod\n",
      "ff-ncaa/NZC.frcmod\n",
      "ff-ncaa/OCY.frcmod\n",
      "ff-ncaa/OMX.frcmod\n",
      "ff-ncaa/ONL.frcmod\n",
      "ff-ncaa/TRO.frcmod\n",
      "ff-ncaa/TY2.frcmod\n",
      "ff-ncaa/TYQ.frcmod\n",
      "ff-ncaa/YCM.frcmod\n",
      "ff-ncaa/YNM.frcmod\n",
      "ff-nucleic-OL15/ff-nucleic-OL15.frcmod\n",
      "ff-ptm/0AF.frcmod\n",
      "ff-ptm/2MR.frcmod\n",
      "ff-ptm/4PQ.frcmod\n",
      "ff-ptm/ALY.frcmod\n",
      "ff-ptm/BTK.frcmod\n",
      "ff-ptm/CGU.frcmod\n",
      "ff-ptm/CSO.frcmod\n",
      "ff-ptm/CSP.frcmod\n",
      "ff-ptm/CSS.frcmod\n",
      "ff-ptm/DA2.frcmod\n",
      "ff-ptm/DAH.frcmod\n",
      "ff-ptm/HY3.frcmod\n",
      "ff-ptm/HYP.frcmod\n",
      "ff-ptm/LYZ.frcmod\n",
      "ff-ptm/M3L.frcmod\n",
      "ff-ptm/MLY.frcmod\n",
      "ff-ptm/ORM.frcmod\n",
      "ff-ptm/P1L.frcmod\n",
      "ff-ptm/PCA.frcmod\n",
      "ff-ptm/PRK.frcmod\n",
      "ff-ptm/PTR.frcmod\n",
      "ff-ptm/SEP.frcmod\n",
      "ff-ptm/TPO.frcmod\n"
     ]
    }
   ],
   "source": [
    "amber.listFiles()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Build and ionize using Amber"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "topos_amber = amber.defaultTopo()\n",
    "frcmods_amber = amber.defaultParam() + [join(datadir, 'benzamidine.frcmod')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-04-30 14:53:04,457 - htmd.builder.amber - INFO - Starting the build.\n",
      "2018-04-30 14:53:20,381 - htmd.builder.amber - INFO - Finished building.\n",
      "2018-04-30 14:53:24,408 - htmd.builder.ionize - INFO - Adding 9 anions + 0 cations for neutralizing and 0 ions for the given salt concentration.\n",
      "2018-04-30 14:53:24,919 - htmd.builder.ionize - INFO - Min distance of ions from molecule: 5A\n",
      "2018-04-30 14:53:24,920 - htmd.builder.ionize - INFO - Min distance between ions: 5A\n",
      "2018-04-30 14:53:24,921 - htmd.builder.ionize - INFO - Placing 9 anions and 0 cations.\n",
      "2018-04-30 14:53:31,947 - htmd.builder.builder - WARNING - Segments ['P0'] contain both protein and non-protein atoms. Please assign separate segments to them or the build procedure might fail.\n",
      "2018-04-30 14:53:32,677 - htmd.builder.amber - INFO - Detecting disulfide bonds.\n",
      "2018-04-30 14:53:32,715 - htmd.builder.builder - INFO - 6 disulfide bonds were added\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bond between A: [serial 349 resid 24 resname CYX chain A segid P0]\n",
      "             B: [serial 569 resid 40 resname CYX chain A segid P0]\n",
      "\n",
      "Bond between A: [serial 1681 resid 115 resname CYX chain A segid P0]\n",
      "             B: [serial 2609 resid 182 resname CYX chain A segid P0]\n",
      "\n",
      "Bond between A: [serial 2144 resid 147 resname CYX chain A segid P0]\n",
      "             B: [serial 2351 resid 161 resname CYX chain A segid P0]\n",
      "\n",
      "Bond between A: [serial 1602 resid 108 resname CYX chain A segid P0]\n",
      "             B: [serial 3002 resid 209 resname CYX chain A segid P0]\n",
      "\n",
      "Bond between A: [serial 2492 resid 172 resname CYX chain A segid P0]\n",
      "             B: [serial 2797 resid 196 resname CYX chain A segid P0]\n",
      "\n",
      "Bond between A: [serial 90 resid 6 resname CYX chain A segid P0]\n",
      "             B: [serial 1986 resid 136 resname CYX chain A segid P0]\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-04-30 14:53:35,928 - htmd.builder.amber - INFO - Starting the build.\n",
      "2018-04-30 14:53:52,280 - htmd.builder.amber - INFO - Finished building.\n",
      "2018-04-30 14:53:55,583 - htmd.molecule.writers - WARNING - Field \"resid\" of PDB overflows. Your data will be truncated to 4 characters.\n"
     ]
    }
   ],
   "source": [
    "bmol_amber = amber.build(smol, topo=topos_amber, param=frcmods_amber, outdir='./build_amber')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Visualize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The built system can be visualized (with waters hidden to be able to visualize the inserted ions):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c848d4df66564809a65ed253fa486a8a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "bmol_charmm.view(sel='not water') # visualize the charmm built system\n",
    "# bmol_amber.view(sel='not water') # uncomment to visualize the amber built system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `bmol_charmm` and `bmol_amber` are `Molecule` objects that contain the built system, but the full contents to run a simulation are located in the `outdir` (`./build` in this case)."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  },
  "widgets": {
   "state": {
    "0a178ab50dc445f7924510f91887f5fc": {
     "views": [
      {
       "cell_index": 15
      }
     ]
    },
    "25c1755bbe564dbab74cd85ab64fd1a9": {
     "views": [
      {
       "cell_index": 25
      }
     ]
    },
    "2b9fe15c349d4639adf691ce9756fdae": {
     "views": [
      {
       "cell_index": 28
      }
     ]
    },
    "574684028db242379374737dca9fbb19": {
     "views": [
      {
       "cell_index": 25
      }
     ]
    },
    "5c49905aacd04d71b82425f44ae9e531": {
     "views": [
      {
       "cell_index": 28
      }
     ]
    },
    "5fec83d14a08443fa1f8a9e24092a958": {
     "views": [
      {
       "cell_index": 34
      }
     ]
    },
    "62600751c0694125afa73b13bfbd1957": {
     "views": [
      {
       "cell_index": 15
      }
     ]
    },
    "d4bbcf3b9c59484aac9988a2f2ff1056": {
     "views": [
      {
       "cell_index": 34
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
